#ifndef GIBBS_HPP
#define GIBBS_HPP
#include "simulate/rtnorm.hpp"
#include "simulate/rtnomNC.hpp"
#include <vector>


class Gibbss
{
	public:
		Gibbss(double s,mat X,double *Y, double xi)
		{
			rt=new rtnormNC();
			_X=X;
			_p=X.n_cols;
			mat v(_p,_p);
			v.fill(s);
			_v=diagmat(v);
			_n=X.n_rows;
			_Y=Y;
			//U=new Distribution::Uniform(1,0,1);
			boost::random::uniform_01<> U;
			runif=new RandomG::Random<boost::random::uniform_01<> >(U);
			_V=inv(inv(_v)+_X.t()*_X);
			_O=_V*_X.t();
			_xi=xi;
		}
		~Gibbss()
		{
			//delete rt;
			//delete U;
			//delete runif;
		}
		mat rz(mat b)
		{
			mat z(_n,1);
			for(int i=0;i<_n;i++)
			{
				double foo=dot(_X.row(i),b);
				if(_Y[i]==0)
				{
					z(i,0)=(*rt)(foo,-pINF,0,_xi);
					//cout << z(i,0) << " y " << _Y[i] << "\n"; 
				}else{
					z(i,0)=(*rt)(foo,0,pINF,_xi);
				}
			}
			return z;
		}
		mat Pi_b(mat z)
		{
			mat m=_V*_X.t()*z;	
			/*cout << V;
			cout << X2;
			cout << nb;*/
			G=new Distribution::Gaussian(_p,m,_V);
			mat res1=(*G).r(1);
			//cout << z.t();
			//cout << beta.t();	
			delete G;	
			return res1;			
		}
		mat Sample(int M)
		{
			mat Gamma(M,_p);
			Gamma.fill(1);
			mat Y=vect2mat<double>(_Y,_n,1);
			mat b=_O*Y;
			mat z=_X*b;
			G=new Distribution::Gaussian(_p,b,_V);
			mat res1=(*G).r(1);
			cout << b;
			delete G;	
			Gamma.row(0)=res1.t();
			for(int i=1;i<M;i++)
			{
				cout << " " << i;
				mat vect=Pi_b(z);
				Gamma.row(i)=vect.t();
				z=rz(Gamma.row(i));
				cout << vect.t();
			}
			return Gamma;

		}
	private:
		rtnormNC *rt;
		Distribution::Gaussian *G;
		mat _v;
		mat _O;
		//Distribution::Uniform *U;
		double *_Y;
		mat _V;
		mat _X;
		int _p;
		int _n;
		RandomG::Random<boost::random::uniform_01<> > *runif;
		mat _Bg;
		double _xi;
		mat _Vg;
};
#endif
